Open source code for the generation of digital reference objects for dynamic contrast-enhanced MRI analysis software validation

Objectives: Dynamic contrast-enhanced MR images can be analyzed through the application of a wide range of kinetic models. This process is prone to variability and a lack of standardization that can affect the measured metrics. There is a need for customized digital reference objects (DROs) for the validation of DCE-MRI software packages that undertake kinetic model analysis. DROs are currently available only for a small subset of the kinetic models commonly applied to DCE-MRI data. This work aimed to address this gap. Methods: Code was written in the MATLAB programming environment to generate customizable DROs. This modular code allows the insertion of a plug-in to describe the kinetic model to be tested. We input our generated DROs into three commercial and open-source analysis packages and assessed the agreement of kinetic model parameter values output with the ‘ground-truth’ values used in the DRO generation. Results: For the five kinetic models tested, the concordance correlation coefficient values were >98%, indicating excellent agreement of the results with ‘ground-truth’. Conclusions: Testing our DROs on three independent software packages produced concordant results, strongly suggesting our DRO generation code is correct. This implies that our DROs can be used to validate other third party software for the kinetic model analysis of DCE-MRI data. Advances in knowledge This work extends published work of others to allow customized generation of test objects for any applied kinetic model and allows the incorporation of B1 mapping into the DRO for application at higher field strengths.


INTRODUCTION
It has previously been documented 1-3 that the complex analysis of dynamic contrast-enhanced MR imaging (DCE-MRI) data using kinetic modeling techniques can involve significant variability, and lack of standardization. This has restricted the translation of DCE-MRI from a research technique into successful application in the clinic. 4 Recently, the Quantitative Imaging Biomarkers Alliance within the RSNA (QIBA) has been established to explore this lack of standardization and other variability, and to establish common protocols and procedures to eliminate it. 5 One source of variability in DCE-MRI kinetic model analysis is the lack of standardized analysis software. There are many software packages available for such analyses, both commercial and open source, and many institutions develop their own bespoke solutions. One author has concluded that software variability, whether due to errors in implementation or variability in design, can be responsible for differences in output kinetic model parameters of up to 74% (within subject coefficient of variation) when presented with identical input data. 6 The availability of suitable input reference data is paramount to the testing and validation of such software packages. In 2009, QIBA made a set of digital reference objects (DROs) available which simulate in silico DCE-MRI input data with known kinetic model parameter output. 7 Each QIBA DRO is delivered as a set of DICOM images that can be used as input into almost any DCE-MRI analysis package.
Although the QIBA DROs are useful, they are only currently available for the Tofts and extended Tofts kinetic models, which limits their wider application. The fact that they have been generated using a third party simulation software package, JSim, 8 which may not be familiar to many readers, makes it difficult to customize or amend these DROs. In addition, they do not account for the full complexity of Objectives: Dynamic contrast-enhanced MR images can be analyzed through the application of a wide range of kinetic models. This process is prone to variability and a lack of standardization that can affect the measured metrics. There is a need for customized digital reference objects (DROs) for the validation of DCE-MRI software packages that undertake kinetic model analysis. DROs are currently available only for a small subset of the kinetic models commonly applied to DCE-MRI data. This work aimed to address this gap. Methods: Code was written in the MATLAB programming environment to generate customizable DROs. This modular code allows the insertion of a plug-in to describe the kinetic model to be tested. We input our generated DROs into three commercial and open-source analysis packages and assessed the agreement of kinetic model parameter values output with the 'ground-truth' values used in the DRO generation. Results: For the five kinetic models tested, the concordance correlation coefficient values were >98%, indicating excellent agreement of the results with 'ground-truth'. Conclusions: Testing our DROs on three independent software packages produced concordant results, strongly suggesting our DRO generation code is correct. This implies that our DROs can be used to validate other third party software for the kinetic model analysis of DCE-MRI data. Advances in knowledge This work extends published work of others to allow customized generation of test objects for any applied kinetic model and allows the incorporation of B 1 mapping into the DRO for application at higher field strengths.
DCE-MRI data acquired in some experimental conditions: for example, the incorporation into the acquisition protocol of a sequence that can generate a B 1 map, is now considered essential when performing DCE-MRI at high field strengths (typically ≥ 3 T). In addition, the customization of the arterial input function (AIF) used in the generation of the QIBA DROs is not easy to undertake.
In summary, there is a current need for publicly available DROs, similar in format to the QIBA DROs, that take account of the whole range of kinetic analysis models commonly in use, including reference B 1 maps as well as baseline T 1 (T 1,0 ) maps (or the equivalent source images from e.g. variable flip angle (VFA) acquisitions) along with the dynamic series itself.
This paper presents open-source MATLAB (Mathworks, Natick, MA) code that is designed to generate configurable DICOM DROs with these attributes. The code is general in purpose and presented in a modular format in order to encourage modification and extension to suit the user's particular requirements. It is primarily designed as a tool for quality control and validation of DCE-MRI analysis software packages, although it could also be used for simulation studies where one or more aspects of the resulting digital phantoms are varied for investigational purposes.

METHODS AND MATERIALS
A generalizable and customizable set of scripts was developed in MATLAB for the generation of DCE-MRI digital reference objects.
The analysis of such DCE-MRI image data via a kinetic model typically proceeds via a number of stages. The first stage is the conversion of dynamic image intensities (gray-scale levels indicative of acquired MR signals) to concentration-time curves, representing the evolution of contrast agent concentration (typically gadolinium based) in each imaged voxel. In a typical DCE-MRI acquisition, this is achieved through correction of the images for imperfect excitation flip angles using B 1 mapping techniques. There then follows computation of the T 1 relaxation time (or equivalently R 1 relaxation rate) evolution over time via the standard SPGR equation 9 that relates the tissue T 1 with the acquired signal intensity. A relaxivity relationship is then applied to yield the required contrast agent concentration time (CACT) curve.
Secondly, a suitable kinetic model is fitted to the CACTs yielded at each voxel and a reference curve, the arterial input function, which is the CACT measured in the artery supplying blood to the tissues in the imaged voxels.
Finally, application of the kinetic model to the CACTs, usually incorporating curve fitting via an optimization algorithm, yields physiologically meaningful parameter values for each imaged voxel.
The algorithm followed to generate DROs reverses the steps above, leading to the construction of synthetic DCE-MRI images from known kinetic model parameter values. This allows subsequent analysis of the synthetic images to compare the output with the known input, and hence to achieve validation of the analysis software used.
MATLAB code was developed to generate customizable DCE-MRI DROs. This code was designed with the following characteristics ( Figure 1): An input domain of kinetic model index values is specified along with parameters defining the layout and composition of the resulting DRO. Also configured are the various parameters associated with the MR acquisition being simulated (e.g., TR, TE, time resolution, duration, etc.) and the synthetic 'tissues' under study (native T 1 values, baseline signal values, relaxivity constants, etc.). In addition, the AIF to be used in the DRO generation is specified as a 'blood-concentration' .
For each voxel in the DRO, the kinetic model under investigation is applied in the 'forwards direction' using the associated known kinetic model indices and the given AIF. The kinetic model is specified via a modular plug-in script containing the residue function specific to the applied model expressed in convolutional form. A contrast-agent concentration-time curve is output from this process, one for each voxel in the tissue to be simulated.
In turn, the CACTs are then converted into MR signal grey-level curves (signal-time curves) by reverse application of the standard relaxivity and spoiled gradient echo (SPGR) equations as implemented by Schabel & Parker. 10 Finally, a 4D dynamic DICOM image set is exported with grey-level evolution in time at each voxel as specified by the signal-time curve.
A T 1 (and the equivalent R 1 ) map is generated for the synthetic tissue volume along with a set of variable flip-angle images which can be used for generation of this T 1 map. Finally, the associated B 1 map, used in calculating the voxel flip-angles for the dynamic series and T 1,0 maps, is also output in DICOM format.
For completeness, the input known value kinetic parameter maps are also output in DICOM format. Suitable DICOM headers are attached to all output images, ensuring that they may be read by most DCE-MRI analysis software packages taking input in DICOM format.
The resulting DRO output volume can be configured to have any size and one or more representative MR image slices.
Commonly used kinetic models, especially in tumor studies, include the Tofts model 11,12 and the extended Tofts model 13 which yield the parameters K trans , v e and, in the latter case, v p . Other models, for example the Patlak model, 14 a tissue uptake model, 15 and a two-compartment exchange model (2CXM), 16 have also been employed where there is an a priori theoretical reason to expect their better representation of the physiology and applied acquisition regime (see [17][18][19] for further details).
Example DROs were generated using our MATLAB code for the Tofts model, the extended Tofts model, the Patlak model, a tissue uptake model and a two-compartment exchange model. The residue function for each was taken from the work of Sourbron & Buckley 17 as implemented in MATLAB in Barnes et al. 20 The arterial input function used was as described in the QIBA DRO documentation 7 and is provided as a '.csv' file with the source code. Application of the model using a residue function requires numerical integration. In our examples, this integration was implemented using a trapezoidal algorithm. To facilitate accurate numerical calculation for DRO's subsequently down-sampled to a low temporal resolution, we allow the user to specify a time resolution for the integration process that may be different from that of the resulting dynamic series. In our examples, both the integration interval and the dynamic series time resolution were set to 0.5 s.
The configuration parameters of each DRO are summarized in Table 1 21 Starting values for the kinetic model optimization process are given in Table 2.
In the case of ROCKETSHIP and MADYM, there is a requirement to input images in NIfTI format. This was achieved using a third party DICOM-to-NIfTI conversion tool ( dcm2niiX. exe) 22 in the case of ROCKETSHIP, and using the supplied conversion tool 'madym_DicomConvert' in the case of MADYM. Sample conversion scripts for the latter case are supplied with the source code.
Agreement between output kinetic model parameter values and the known true values was displayed in the form of concordance plots, and expressed numerically in terms of the concordance correlation coefficient (CCC) 23 with 95% confidence intervals. Examples of non-perfect concordance were investigated by running the analysis software packages with differing starting values for each of the modeled parameters.
Additional DROs were also generated and tested for the eTM and 2CXM models using different dynamic series time resolutions and integration intervals, to investigate any resolution dependent effects.

RESULTS
The essential features of one example DRO generated by our code are shown in Figure 1.
Sample results are shown in Figures 2 and 3 as concordance plots (drawn using a custom MATLAB script) and expressing the agreement between the kinetic parameter maps output by each analysis platform, and true parameter values (i.e., those used to generate the DROs in the first place). Table 3 shows CCC values for all parameters for all models in the three analysis packages tested. The CCC values are >99% in all but one case which indicates almost perfect agreement 24 with the known true values in every model and software package. Note that for the eTM, voxels with a true K trans of 0.0 are excluded from the plotted data in Figure 2. Convergence of the model was typically poorer for this subset of voxels representing a largely non-physical case.
We can conclude from this that, in all likelihood, each software package is correct and that the DROs generated by our code are correct.
The CACT curve and the fitted curve were plotted for each of the data-points not lying on the concordance line in Figures 2 and 3. It was found that perfect, or in a minority of cases near-perfect, concordance could be restored by changing the initial values of the fitted parameters used in the optimization process. This would indicate that local minima found in optimization were the primary cause of non-perfect concordance of parameter values. Further detail on this can be found in the Supplementary Material 1. In Figures 2 and 3 the number ('n') of distinct parameter combinations analyzed is fewer in ROCKETSHIP than in MADYM or MIStar since the former cannot include a B 1 map in the analysis process.

DISCUSSION
Open source MATLAB code has been presented here for the generation of customizable and generalized digital reference objects for the kinetic model analysis of DCE-MRI data. This code is offered for the production of validation in silico phantoms, although this use could be extended to more detailed simulation studies. Five sample DROs were generated and shown to be correct in the sense that they yielded results in excellent agreement with the known input values when analyzed by three separate third-party analysis software packages.   There are limitations to our DRO study and code base. We are not aware of any true gold standard in DCE-MRI kinetic model analysis. Although physical perfusion phantoms have been developed 32,33 to validate CACT curves calculated from DCE-MRI experiments, these do not directly validate any known kinetic parameter associated with the phantom (e.g., K trans ).
The use of kinetic model code in generating a DRO inevitably leads to the possibility of variability in the DRO images according to the algorithms employed. Since the coding of the kinetic model algorithms in our DRO code follows the ROCKETSHIP implementation, they might be expected to yield good validation results when analyzed using the ROCKETSHIP software. This is less expected in the case of the MIStar and MADYM software packages where details of the kinetic model implementation are either unpublished (MIStar) or independently implemented (in C ++ in MADYM) and hence not studied by the authors of this paper.
It was found that the degree of agreement of results with the known true values was to some extent dependent on the starting values of the kinetic model parameters in the associated  optimization process. This is as expected in general in optimization processes and was especially apparent for the 2CXM case with four free parameters.
With the proviso that the DROs were generated using a dynamic series (and AIF) time resolution of less than 5 s, they were found to yield excellent concordance with known values when analyzed by three independent software packages. The dependence of analysis results on the time-resolution of the dynamic series (and especially of the AIF) is to be expected and was more pronounced in the 2CXM than the eTM, no doubt again due to the inherent complexities of fitting a greater number of free parameters. Interestingly, the ROCKETSHIP analysis showed no degradation of CCC values with dynamic series time resolution when the integrals were calculated with a matched time interval. We can hypothesize that this is because the ROCKETSHIP algorithm was most nearly matched to the DRO generation algorithm in these circumstances (i.e., trapezoidal integration at the time resolution of the dynamic series). MADYM showed pronounced failure of CCC in the same (poor) DRO generation conditions indicating that its analysis algorithm may more correctly interpolate the CACT curves before integration. When the DRO was generated more correctly using a small integration time interval, MADYM performed with much greater concordance. The user of our software should therefore select time resolution of generated DROs with care: a dynamic series/AIF time resolution of ≤2 s and an integral time resolution of ≤0.5 s is recommended for the most consistent results.
It can be noted that the synthetic DCE-MRI images produced by our code are free from noise or artefacts and can be relied upon to yield kinetic parameter indices without complications arising from typical imperfections present in actual in vivo imaging. This is both an advantage, in that any error in analysis is therefore not obscured or confounded by errors in acquisition, and a disadvantage in that the scope for realistic simulation studies is diminished. A rudimentary capability for adding noise to the generated DRO images is provided in our software although realistic simulation is not the primary aim of this study. In generating the DROs we have applied the kinetic model to the CACT curves rather than the signal intensity time curves themselves. This could introduce unwelcome non-linearities in error propagation in signals that include noise. Therefore, use of these DROs for simulation studies would not especially be recommended; comprehensive simulation studies would be coded much more efficiently through other means.
The user may need to edit or augment the supplied DICOM header information for any particular analysis application. This can be accomplished easily in the MATLAB programming environment or indeed using many commonly available DICOM viewers.
Conversion of our DROs to NIfTI format 34 is also easily accomplished using third party DICOM-to-NIfTI converters. The reverse conversion (NIfTI-to-DICOM) is seldom so straightforward, which is the reason we chose to export our DROs as DICOM objects.

CONCLUSIONS
We have presented here a substantially validated tool for quality control and evaluation of DCE-MRI kinetic model analysis packages. Our code to generate DROs is in a widely known format (MATLAB), and is easily generalizable and customizable. We also make available in DICOM format some example DROs constructed for common kinetic models in clinically realistic scenarios.